Statistical mechanics of double-helical polymers 
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We introduce a simple geometric model for a double-stranded and double-helical polymer. We study the 
statistical mechanics of such polymers using both analytical techniques and simulation. Our model has a single 
energy-scale which determines both the bending and twisting rigidity of the polymer. The helix melts at a partic- 
ular temperature Tc below which the chain has a helical structure and above which this structure is disordered. 
Under extension we find that for small forces, the behaviour is very similar to worm-like chain behaviour but 
becomes very different at higher forces. 
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Recent developments in single molecule manipulation tech- 
niques has led to the detailed study of the mechanical prop- 
erties of DNA and other biomolecules as well as their re- 
sponse to applied fields. The model most used in the study 
of the large-scale properties of biopolymers is that of the 
worm-like chain [jl|] in which the polymer flexibility (struc- 
ture) is determined by a single length, the persistence length 
Lp which measures the tangent-tangent correlations. For ex- 
ample, DNA has a persistence length Lp w 50nm. Such 
coarse-grained models are needed to understand the statistical 
behaviour of long chains with a large number of monomers. 
They are complementary to chemically specific models which 
describe accurately the local effects of external fields but gen- 
erally cannot deal with long chains Generalisations of the 
worm-like chain to introduce twist degrees of freedom have 
also been proposed [||] . Whilst the worm-like chain model and 
its generalisations give a good account of DNA under small 
fields (perturbations), it fails when these perturbations become 
large, e.g. under tensional forces above 65pN Such sit- 
uations are biologically relevant in for example in DNA repli- 
cation and repair 

DNA is double-stranded and helical. We introduce a ge- 
ometric model which includes three important ingredients of 
DNA in B-form. First, the double-stranded nature of DNA 
given by the two phosphate backbones; second, the hydrogen 
bonds that keep the two strands together and third, the interac- 
tions between the bases that bring about the twisted structure. 
Similar models could be used to study other double-helical 
polymers Uke F-actin. The model is a generalisation of the 
double-stranded semiflexible (ribbon) polymer introduced by 
Liverpool, Golestanian and Kremer (LGK) [^] which takes 
into account the first two aspects but ignores the third. Other 
ribbon models have also been proposed by a number of au- 
thors iQ]. It will be our conclusion that it is exactly this third 
property, the base-stacking interaction, which can account for 
many of the elastic properties of DNA. Unlike Z. Haijun et 
al who proposed another generalisation of the LGK model 
with base-stacking interactions, we find that the rigidity of 
the individual strands is irrelevant for the effective persistence 
length, Lp. This is consistent with the large difference be- 



tween the persistence length of single and double-stranded 
DNA measured experimentally. Our model can easily be ex- 
tended to include excluded volume and electrostatic interac- 
tions which for the moment have been ignored . We do not 
attempt to make quantitative comparison with experiments but 
point out qualitative differences between simpler models. Our 
aim is rather to suggest a minimal model required to under- 
stand experimental results. 




FIG. 1. (a) The schematic double-helical polymer with 
R2 (s) — Ri (s) + ab(s). The mid section (solid) is in the minimum 
energy configuration, (b) The discretisation of the double-helix used 
in the simulation. 

The polymer is embedded in three dimensional space. 
As for the ribbon, we see qualitatively different be- 
haviour for the high and low temperature regimes with 
a melting of the helix at a particular temperature Tc ~ 



B/kBd where P and a are respectively 



'4\2 r 5+(2P/7ra)n ^/^ 
\9) [l + (2P/7ra)^^ 

the pitch and diameter of the helix at T = and B the 'stiff- 
ness' of the base-stacking interaction. Below Tc, there is a 
helical structure of the chain whilst above the helical struc- 
ture vanishes. Unlike the ribbon, both the bend and twist rigid- 
ity of the double-helical polymer are determined by the base- 
stacking interaction and the stiffness of the individual strands 
is irrelevant. We also study the response of the double-helical 
polymer to an extensional force and find that for small forces, 
the mechanical response is very similar to that of a worm-like 
chain but becomes significantly different at higher forces. 
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The system is composed of two semi-flexible chains mod- 
elUng the sugar phosphate backbones, each with rigidity A, 
whose embeddings in 3-dimensional space are defined by 
Ri(s) and R2(s). Because of the covalent bonds (Ecb ^ 
lOOksT at room temperature) between the ahernating phos- 
phate and sugar molecules, the backbones are effectively inex- 
tensible at the forces and temperatures we consider Similarly, 
the H-bonds {EHb ^ lOfcsT) between the two backbones are 
assumed to keep the distance between the chains constant. 

The separation of the strands is fixed and given by a, and the 
H-bonds give the constraint R2(s) = Ri(s) + ab(s) where 
b(s) is a unit vector which we call the bond-orientation vec- 
tor (see Fig. 0). The bond-orientation vector is perpendicular 
to both strands. The infinitesimal distance in 3-dimensional 
space between neighbouring points on a single strand is a 
function of the local curvature , 9^Rj(s) of the strands, where 
dxA{x) = dA/dx. The base-stacking interaction can be 
modelled by a potential Vi(9gRi,b), where the subscript 
i G {1,2} refers to strand 1 or 2, which is a function of the 
curvatures of the chains and the bond-vector only, with a min- 
imum given by a symmetric double-helix. The Hamiltonian 
of the system is given by 



n = y"ds{Vi(92Ri,b) + y2(9'R2,b) 



+A(|5fRip + |a2R,|2)} 



(1) 



The simplest example of such a potential which we will use 
for the remainder of this letter is the quadratic Vi (9^Ri , b) = 
Bid^Ri-Hh)^ and ¥2(8^112, b) = ^(a^Ra + iJb)^. This 
is a potential of stijfness B with a minimum when equivalent 
points on the two strands have equal and opposite curvatures, 
±iJ. More complicated potentials with several minima are 
also possible (see later). We note that we have not included 
the asymmetry of the DNA helix but a similar but more com- 
plicated model could be defined for an asymmetric double- 
helix [^. We implement the model by introducing the "mid- 
curve" r(s): Ri(s) = r(s) + f b, R2(s) = r(s) - f b. 
In terms of the tangent to the mid-curve t(s) — dgV and the 
bond-director b the Hamiltonian of the system can now be 
written as 



[dstis) + ^d^s) - Hh 



(2) 



subject to the exact (local) constraints 



(t ± ^ dshr 



1, 



1, 



(t±-a3b).b = o. (3) 



corresponding to the chain inextensibility, fixed distance be- 
tween strands and definition of the bond vector respectively. 
We make the assumption that B ^ A since the bending rigid- 
ity of single-stranded DNA is very smafl, (A/kBT < lOA), 
compared to the double-stranded form. We can define a typi- 
cal length £ = B/ksT, which together with the strand sepa- 
ration a and equihbrium radius of curvature iJ^^/^ form the 



three relevant lengths of the problem. This completes the for- 
mulation of the model. 

T=0 behaviour We can calculate the ground state of 
this system from the condition H ^ which gives, 

Ri(s) = f cos (Gs) ei + f sin (6 s) ea + s-y/l - and 
R2(s) = -f cos (6 s) ei - § sin (9 s) 62 + s^Jl -^83 
where O = \ —■ This is a double helix with a pitch: 
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and the bond-director field, b(s) = 
— cos (9 s) ei — sin (9 s) 62 where {ci} are an orthonormal 
set of unit vectors. 

Finite T behaviour Inspite of the simplification of the 
model from using constraints, an exact analytic expression 
for the partition function of the model is still not avail- 
able. For our analytic calculations, we use a very useful 
approximation [|o[ which describes correctly the equilib- 
rium behaviour of correlation functions and even distribu- 
tion functions, if the end effects are taken into account care- 
fully [|ri]|. This is a mean field approach that relaxes the 
local constraints to global ones [ p^HTl] ]. It may be thought 
of as a self-consistent theory and corresponds to the saddle- 
point evaluation of path integrals over the lagrange multi- 
pliers [|l2|] . This approach fails for some dynamical proper- 
ties [|13[] which are outside the scope of this letter. To im- 
plement this, we add a variational term to our Hamiltonian 
^ = /ds[|(t-f W + |(t + f W + ||^b2 

+f (t - § dsh) • b + f (t + f dsh) ■ b] , where b, c, and e 
are dimensionless constant Lagrange multipliers . 

Two useful geometric parameters u — Ha/2 and v — 
4(£/a)^ can be used to characterise the double-helical poly- 
mer. Note that v is proportional to T^^ and can be viewed as 
a measure of temperature. We then determine the constants 
self consistently by demanding the constraints of Eq.(|^) to 
hold on average, where the thermal average is calculated by 
using the total Hamiltonian H + Hm- The self-consistency 
lead to the following set of equations for the constants b 



and c: (c + u^v'^ 



UV) + Vc + M^W^ I = 



9v^ 



32 ' 4%/26 ^ 

0. The above equations, which 



, V ^ I — k and e 

3 V ^ 3 

are nonlinear and difficult to solve exactly, determine the be- 
haviour of b and c functions of u and v. One can solve the 
equations analytically in some limiting cases corresponding 
to T and T -> cx). For u ^ cx) (T ^ 0) and 7^ we 



find, b 
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UV and for v 



(T — > 00) and H ^ 0,we obtain b ^ |, c 
A full solution requires a numerical treatment and shows sim- 
ple monotonic behaviour for both c and b jl5| ] . 

We can then calculate the correlation functions. For the 
tangent-tangent correlation one obtains 



(t(.).t(O)) 



-\/2fc| 



T->0 



4V2& 



(1 



(4) 



whereas for the bond-director field one obtains 



2 



(b(.).b(O)) 



3z e 
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T-+0 



{uv -b) + VD£^ ^ {uv ~ b) ~ VdP _ 



e »" f cos 



(5) 



We have defined 



discriminant D{u,v) 



(c + u^v^). The tangent-tangent correla- 



tion (Eq.(^) of the double-helical polymers is similar to 
that of a worm-like chain but with a length rescaling fac- 
tor, and decays exponentially at all temperatures. Eq.(|]) on 
the other hand, indicates a change of behaviour at L* = 
i.e. (6 — Mw)2 = c + u^w^ for the bond-director correlation. 
The correlation is over-damped for {b — uv)"^ > c + v^v^ 
(high temperatures), while it is under-damped (oscillatory) 
for (6 — uvY < c + u^v^ (low temperatures). The interesting 

point (6 — uv)"^ = c + v^v^ happens for Wc — 9 (^ 4(1^4^-) ^ 
that leads to the value for Tc quoted above. We also find 
a divergence in the specific heat, Cv at Tc- We emphasise 
that it is not a thermodynamic phase transition in the sense of 
long-range ordering and broken symmetry. It is a cross-over 
that appears due to competing effects, and the transition is 
from a state with some short-range order to a state with a 
different short-range order The cross-over (transition) point 
corresponds to a type of 'Lifshitz point' for a 1-d system JTil]. 
At high temperatures, in addition to the effects we describe we 
also expect the H-bonds between the bases to break, leading 
to denaturation which we cannot treat in this model. In the 
low T regime, we can read off the bend (tangent) persistence 
length Ltp — — and twist (bond) persistence 

length, isp = l^i/a and helical frequency (2iJ/a)^/'^. Both 
the bend and twist rigidity of the double-helical polymer are 
given by the same energy scale B = ksTi and have little to 
do with the stiffness of the individual strands. 

Force-extension: unwinding the helix Under a stretching 
force, which without loss of generality we orient along the 
z-axis F = Fi, the Hamiltonian becomes T-Lf — Ti. — F ■ 
J dsdsT. Now (tz) — (t • z) ^0 and we must use con- 
nected correlation functions to calculate the self-consistent 
equations. Defining / = Fa/8 we obtain the equations 



4V26 
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(c + u\^) [{b - uv) + Vc + u^v^] = 



(6) 



We can solve eqns. (^ for b{u, v, /) and c{u, v, /) and the 
extension per unit length z{u,v, f ) is simply: 



z = -^^ {t{s))ds = ^^/v 



(7) 



We can therefore obtain the force-extension relation for the 
double-helical polymer Eq. ^ has been solved numerically 
and in Fig. ^ we have plotted z against /. In the inset we 
compare the force-extension curve to that of the WLC and 
find they are identical for small forces. We note the interest- 
ing fact that the mean field model can also be used to calculate 



the force-extension behaviour of the WLC. It has the simple 
compact form fLp — -^i^^r^rji- We find that this simple ex- 
pression has an error of less than 4 % compared to the 'exact' 
value of Marko and Siggia [pl|. 




FIG. 2. Force-extension relation for the double-helical polymer 
(DHP). This is a log plot with u = 0.77 and w = 115000, typical 
values for B-DNA which corresponds to a low temperature T < Tc 
(Note / = Fa/8).The dashed line shows the T = Q contour length 
^1 - a2eV4. Inset: comparison with worm-like chain at small 
forces. 

Simulations We also performed Monte Carlo simulations 
of the double-helical polymers. We used a square lattice of 
size A to discretise the 'helical ribbon' and to implement all 
the constraints (see Fig. |l]). The discretised Hamiltonian is 
given by 



Af-l 




-(-i)-'"iib„.«+i-t^„; 



(8) 



where n is the monomer label (position along the chain) and 
j £ {1,2} the strand label (see Fig. ^. Since the Hamiltonian 
is local, we grew chains with a local algorithm. In our simu- 
lation the number of monomers N was taken to be 1000 for 
each chain, A = 1 and a — 3. Typical conformations at high 
and low temperatures are shown in Fig. 0. 

Typical bond correlation functions are plotted in Fig.|| 
From the bond and tangent correlation functions, we can ob- 
tain the bond and tangent persistence lengths. We calculate 
the bond correlation length as Lbp = Cbu£ and the tangent 
correlation length as Ltp = Ct£ with the constants given 
Cb = 1.09 ± 0.09 and Ct = 2.46 ± 0.73 respectively. A 
comparison with the mean field results shows that the results 
are the same to within a constant. This difference can be con- 
sidered as due to finite corrections to the mean field solution 
due to higher order terms in the 1 /d expansion. 
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FIG. 3. The (b(s) ■ b(0)) correlation function measured in the 
simulations for temperatures £ = 400, 50, 3 and u — 0.6 corre- 
sponding to r < Tc and T > Tc- The averages were done over 200 
statistically independent samples. 



forces the extension changes abruptly and increases steeply as 
the double helical structure is unwound. This qualitatively 
agrees with the experimental results on over-stretched DNA 
which nevertheless have a sharper transition and a flat- 
ter plateau. The graph is plotted on a log scale, so as to make 
clearer the complex behaviour of the force-extension curve. 
There are a two obvious ways in which a sharper transition 
could be obtained. First, we have used a simple quadratic 
potential for the base-stacking interaction corresponding only 
to an expansion about the minimum energy conformation. A 
more realistic potential with a cut-off and/or several minima 
could change the sharpness of and number of transitions. In 
addition one could imagine as in mechanical unzipping exper- 
iments on DNA, that there could be cooperative effects due to 
the effects of sequence disorder 

We have benefited from discussions with R. Everaers, R. 
Golestanian and K. Kremer. 



In conclusion, we have studied the properties of a well- 
defined model of a double-helical, double-stranded semi- 
flexible polymer using a mean-field analytical approach as 
well as extensive MC simulations. We have shown novel non- 
trivial differences between the high, low and zero temperature 
behaviour A detailed comparison with the simulations and 
the mean-field solution shows qualitative and almost quantita- 
tive agreement. 




FIG. 4. Typical conformations at (a) low, (b) intermediate and (b) 
high temperatures 

The fine structure studied in this model gives a different be- 
haviour compared to a WLC. We see this most clearly when 
the double-helical polymer is subject to an external pulling 
force. It shows two distinct regimes. For low forces the dy- 
namical response is simular to a WLC, with the most relevant 
energetic contribution given by the bending rigidity. At higher 
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